# install.packages('sf')
library(sf)6 Trabajando con datos espaciales en R
Geoestadística i Tècniques d’Observació Global
6.1 Introducción
El objetivo de este \(lab\) es familiarizarse con el uso y representación de información espacial en \(R\). Además de las funciones estadísticas que ya hemos visto, existen en \(R\) miles de paquetes, cada uno con un propósito específico. Varios de ellos permiten interactuar con datos espaciales, tanto en formato vectorial como raster. Los principales paquetes “espaciales” son:
- rgdal
- sf
- terra
- stars
- tmap
Estas librerías (paquetes) permiten leer y gestionar objetos espaciales, ya sea en formato vectorial o ficheros raster como ASCII o cualquier otro formato soportado por GDAL. Estos paquetes funcionan como un GIS, permitiéndonos realizar la mayoría de geoprocesos (intersecciones, uniones, etc.), así como procesos geoestadísticos avanzados y crear mapas.
Durante años, tres de los paquetes más populares para trabajar con datos espaciales en R han sido rgdal, raster y sp, que se convirtieron prácticamente en un estandard. Sin embargo, por diversos problemas, en 2022 se decidió deprecarlos, es decir, dejar de darles soporte, y se recomendó sustituirlos por sus “equivalentes”: terray sf. Desde octubre de 2023, aún pueden usarse las funciones de sp y raster, pero rgdal ha sido retirado, por lo que nos centraremos en conocer las funcionalidades de terray sf.
Obviamente, no podremos profundizar en toda la potencialidad de estos dos paquetes, que cada vez más logran sustituir perfectamente a un GIS de escritorio, pero sí veremos algunas de sus opciones más básicas.
6.2 Trabajando con vectores
6.2.1 Cargar información vectorial
Sabemos que los archivos shapefile de ESRI, que son casi un estándar en el mundo del SIG, son en realidad un conjunto de varios ficheros. Pero para cargarlos en \(R\), al igual que hacemos en ArcMap, basta con cargar el que tiene extensión .shp, ya que el resto van asociados a él. Sin embargo, si queremos copiar o cortar nuestros “shapes” en otra carpeta debemos tener cuidado de copiar todos los ficheros.
Veamos un ejemplo, importando las capas estaciones_meteo.shp y provincias.shp, que se encuentran en la carpeta datos/shapes, dentro de la sección “Recursos” del campus virtual de la asignatura. Estas capas contiene la información espacial de parte de las estaciones meteorológicas de las que extrajimos la información para ajustar la regresión lineal, así como los límites de todas las provincias de España.
Para leerlos en R podemos usar la función st_read() del paquete sf, por lo que antes que nada, instalamos y cargamos el paquete sf:
Y posteriormente cargamos los ficheros, especificando la ruta al archivo .shp y asignándoles un nombre:
estaciones <- st_read('data/meteo/meteo_espacial/estaciones_meteo.shp') Reading layer `estaciones_meteo' from data source
`C:\Users\Usuari\OneDrive - udl.cat\Teaching\EFCN_Geoestadística\Labs_Geoestadistica_Forestal\data\meteo\meteo_espacial\estaciones_meteo.shp'
using driver `ESRI Shapefile'
Simple feature collection with 90 features and 13 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 574966 ymin: 4514365 xmax: 818327 ymax: 4635648
Projected CRS: ETRS89 / UTM zone 30N
Vemos que, al leer los ficheros, la consola nos informa sobre el nuevo objeto que estamos creando. Identifica que el fichero de origen es un shapefile de ESRI, y nos informa de que se convertirá en un objeto “simple feature” con 90 observaciones de tipo “punto” y 13 campos. También nos informa sobre el bounding box del objeto espacial (el rectángulo que lo contiene) y su sistema de referencia de coordenadas o CRS.
Un objeto simple feature es el método que ha elegido sf para lidiar con información espacial. Los objetos espaciales almacenan tanto información relativa a las características espaciales (sistema de referencia, extensión, coordenadas de los objetos…) como a los atributos de cada uno de los objetos (en este caso, la información meteorológica). La ventaja de sf es que lo hace de manera muy lógica y sencilla. Vemos en el panel de Environment que estaciones aparece como un data frame normal. Y es que de hecho es un data frame normal, pero ahora es a la vez más cosas:
class(estaciones)[1] "sf" "data.frame"
La diferencia es que contiene, además de los atributos que veríamos en ArcMap, una columna llamada geometry con las coordenadas de cada una de las observaciones. La ventaja de esto es que el objeto espacial es un data frame, y por lo tanto podemos procesarlo con normalidad, filtrando, creando nuevas variables o modificando valores, igual que hacemos con los data frames. Al cargar estaciones nos informaba de que se traba de un objeto espacial de puntos. Por tanto cada punto contendrá las coordenadas x e y de ese punto según el CRS del objeto espacial:
estacionesSimple feature collection with 90 features and 13 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 574966 ymin: 4514365 xmax: 818327 ymax: 4635648
Projected CRS: ETRS89 / UTM zone 30N
First 10 features:
OBJECTID FID_weathe INDICATIVO AÃ_O MES NOMBRE ALTITUD
1 1 12 9382E 2012 6 PANCRUDO (D.G.A.) 1285
2 2 85 9987N 2012 6 DELTEBRE (PARC NATURAL) 4
3 3 11 9377E 2012 6 EL PEDREGAL 985
4 4 50 9567U 2012 6 EJULVE (D.G.A.) 1095
5 5 39 9530V 2012 6 UTRILLAS (D.G.A.) 960
6 6 0 3013 2012 6 MOLINA DE ARAGON 1063
7 7 70 9939A 2012 6 FUENTESPALDA (DGA) 720
8 8 40 9531X 2012 6 MONTALBAN 'AUTOMATICA' 885
9 9 84 9981A 2012 6 TORTOSA (OBSER. DEL EBRO) 48
10 10 89 9999 2012 6 ODON 1110
T_MAX_abs T_MIN_abs TMed_MAX TMed_MIN Tmed_MES Provincia
1 29.5 3.0 22.0 9.6 15.80 Teruel
2 32.0 10.0 27.5 17.4 22.45 Tarragona
3 31.0 3.5 21.5 10.2 15.85 Guadalajara
4 30.0 5.0 21.9 11.5 16.70 Teruel
5 32.5 3.5 23.7 10.8 17.25 Teruel
6 32.8 3.2 23.8 8.4 16.10 Guadalajara
7 32.0 5.0 24.4 13.2 18.80 Teruel
8 31.2 5.4 23.0 11.0 17.00 Teruel
9 34.4 13.0 29.3 17.1 23.20 Tarragona
10 32.0 4.0 22.4 10.2 16.30 Teruel
geometry
1 POINT (666121 4514365)
2 POINT (814489 4514857)
3 POINT (620767 4515276)
4 POINT (705449 4516865)
5 POINT (681326 4520030)
6 POINT (593975 4522166)
7 POINT (758868 4522277)
8 POINT (686217 4522281)
9 POINT (794466 4524785)
10 POINT (620133 4526863)
Vamos a cargar ahora la capa provincias:
provincias <- st_read('data/meteo/meteo_espacial/provincias_spain.shp') Reading layer `provincias_spain' from data source
`C:\Users\Usuari\OneDrive - udl.cat\Teaching\EFCN_Geoestadística\Labs_Geoestadistica_Forestal\data\meteo\meteo_espacial\provincias_spain.shp'
using driver `ESRI Shapefile'
Simple feature collection with 51 features and 12 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -14129.47 ymin: 3892590 xmax: 1126923 ymax: 4859517
Projected CRS: ETRS89 / UTM zone 30N
El proceso es idéntico, pero ahora nos dice que el objeto creado contiene una geometría tipo multipolygon, es decir, que es un fichero vectorial de polígonos. Ahora, por lo tanto, el campo geometry contendrá, para cada observación, el conjunto de coordenadas que define los vértices del polígono:
provinciasSimple feature collection with 51 features and 12 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -14129.47 ymin: 3892590 xmax: 1126923 ymax: 4859517
Projected CRS: ETRS89 / UTM zone 30N
First 10 features:
OBJECTID INSPIREID COUNTRY
1 1 ES.IGN.BDDAE.34160100000 ES
2 2 ES.IGN.BDDAE.34080200000 ES
3 3 ES.IGN.BDDAE.34100300000 ES
4 4 ES.IGN.BDDAE.34010400000 ES
5 5 ES.IGN.BDDAE.34070500000 ES
6 6 ES.IGN.BDDAE.34110600000 ES
7 7 ES.IGN.BDDAE.34040700000 ES
8 8 ES.IGN.BDDAE.34090800000 ES
9 9 ES.IGN.BDDAE.34070900000 ES
10 10 ES.IGN.BDDAE.34111000000 ES
NATLEV
1 https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/3rdOrder
2 https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/3rdOrder
3 https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/3rdOrder
4 https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/3rdOrder
5 https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/3rdOrder
6 https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/3rdOrder
7 https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/3rdOrder
8 https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/3rdOrder
9 https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/3rdOrder
10 https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/3rdOrder
NATLEVNAME NATCODE NAMEUNIT CODNUT1 CODNUT2 CODNUT3 Shape_Leng
1 Provincia 34160100000 Araba/Álava ES2 ES21 <NA> 639475.2
2 Provincia 34080200000 Albacete ES4 ES42 <NA> 770993.4
3 Provincia 34100300000 Alacant/Alicante ES5 ES52 <NA> 583497.3
4 Provincia 34010400000 Almería ES6 ES61 <NA> 632690.0
5 Provincia 34070500000 Ávila ES4 ES41 <NA> 630836.7
6 Provincia 34110600000 Badajoz ES4 ES43 <NA> 1136844.0
7 Provincia 34040700000 Illes Balears ES5 ES53 <NA> 1459098.3
8 Provincia 34090800000 Barcelona ES5 ES51 <NA> 850790.3
9 Provincia 34070900000 Burgos ES4 ES41 <NA> 1138046.5
10 Provincia 34111000000 Cáceres ES4 ES43 <NA> 965487.3
Shape_Area geometry
1 3035458537 MULTIPOLYGON (((519021.3 47...
2 14920858305 MULTIPOLYGON (((539277 4215...
3 5820590272 MULTIPOLYGON (((697663 4195...
4 8767315137 MULTIPOLYGON (((496730.8 39...
5 8047820788 MULTIPOLYGON (((292983.8 44...
6 21793269761 MULTIPOLYGON (((165751.6 42...
7 5018967692 MULTIPOLYGON (((868245.7 43...
8 7762160343 MULTIPOLYGON (((898327 4573...
9 14279022189 MULTIPOLYGON (((419298.4 46...
10 19885694014 MULTIPOLYGON (((172175.6 43...
provincias$geometry[1]Geometry set for 1 feature
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 476582.6 ymin: 4702304 xmax: 562700.7 ymax: 4784934
Projected CRS: ETRS89 / UTM zone 30N
MULTIPOLYGON (((519021.3 4717987, 518976.3 4717...
6.2.2 Representando gráficamente información vectorial
Evidentemente, cuando trabajamos con datos espaciales una de las acciones más interesantes es representarlos gráficamente. Para ello basta con usar la función plot() y el nombre del objeto espacial
plot(estaciones)Warning: plotting the first 9 out of 13 attributes; use max.plot = 13 to plot
all

Por defecto nos ploteará todos los campos. Si queremos visualizar específicamente alguno de ellos debemos indicarlo expresamente, poniendo el nombre del campo entre corchetes []:
plot(estaciones["Tmed_MES"])
\(R\) mapea los puntos según su ubicación, y les asigna un color en función de los valores de la variable elegida, pero podemos editar la visualización igual que hacemos con los \(scatterplot\) o los \(plots\) normales.
Por ejemplo, podemos cambiar el tipo de símbolo con pch, el color con col o el tamaño con cex.
plot(estaciones["Tmed_MES"], pch = 19)
plot(estaciones["Tmed_MES"], pch = 21, col='black', bg= 'red')
En caso de que sólo queremos plotear la geometría del objeto (sin incluir ningún atributo) debemos indicarlo con la función st_geometry()
plot(st_geometry(provincias))
plot(st_geometry(provincias), col = "dark red")
Podemos también combinar dos capas diferentes mediante el comando add = TRUE:
plot(st_geometry(provincias) )
plot(estaciones["Tmed_MES"], pch = 19, col = "red", cex = 0.5, add = TRUE) 
6.2.2.1 Visualizando capas vectoriales con tmap
Un paquete interesante para la visualización de información espacial es tmap. Este paquete funciona de manera diferente, ya que le debemos indicar las diferentes capas a visualizar uniendo las órdenes con el símbolo +. Lo primero es indicar la capa que queremos visualizar mediante el comando tm_shape() y luego especificamos si queremos ver puntos (tm_dots()), los bordes de los polígonos (tm_borders()), o los polígonos con relleno (tm_polygons(), entre otras muchas opciones). De cada uno de los comandos tm_* puede personalizarse el color, forma, tamaño, etc. con valores fijos o en función de alguna columna.
Probemos a visualizar las estaciones:
library(tmap)Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
remotes::install_github('r-tmap/tmap')
tm_shape(estaciones) +
tm_dots()
Lógicamente, podemos cambiar el tamaño de los puntos, su forma, o los colores, haciéndolos depender de la variable Tmed_MES e incluso definir la paleta:
tm_shape(estaciones) +
tm_dots(size = 1, shape = 20, col = "Tmed_MES", palette = "viridis")
Con tmaps podemos visualizar tantas capas como queramos, sólo necesitamos volver a utilizar la función tm_shape() con el nombre de la nueva capa. Eso sí, la extensión del mapa vendrá dada por la primera capa que llamemos. Probemos a visualizar las estaciones sobre las provincias del nordeste:
tm_shape(provincias) +
tm_polygons(col = "lightgrey") +
tm_text("NAMEUNIT") +
tm_shape(estaciones) +
tm_dots(col = "Tmed_MES", palette = "viridis", size = 0.75)
La extensión espacial del mapa vendrá determinada por la primera capa o shape que carguemos. Si cargamos primero las estaciones:
tm_shape(estaciones) +
tm_dots(col = "Tmed_MES", palette = "viridis", size = 0.75) +
tm_shape(provincias) +
tm_borders() +
tm_text("NAMEUNIT") 
Las opciones generales del mapa se pueden personalizar mediante tm_layout():
tm_shape(estaciones) +
tm_dots(col = "Tmed_MES", palette = "viridis", size = 0.75,title = "T media") +
tm_shape(provincias) +
tm_borders() +
tm_text("NAMEUNIT") +
tm_layout(legend.outside = F, bg.color = "steelblue",title = "Estaciones meteorológicas del Nordeste", title.size = 4)
Incluso podemos definir un fondo basado en un proveedor como Google Maps o OpenStreetMap, con la función tm_basemap() y hacer que el mapa sea interactivo definiendo antes tmap_mode("view"): (una lista de las opciones se puede encontrar aquí)
tmap_mode("view")tmap mode set to interactive viewing
tm_shape(estaciones) +
tm_dots(col = "Tmed_MES", palette = "viridis") +
tm_shape(provincias) +
tm_borders() +
tm_text("NAMEUNIT") +
tm_basemap("OpenStreetMap.Mapnik")